We bouwen een populatie met een set aan variabelen waarvan het gemiddelde en spreiding bekend is. Hieruit gaan we dadelijk steekproeven (samples) trekken en onderzoeken welk effect de omvang van de steekproef heeft op de nauwkeurigheid waarmee we een uitspraak kunnen doen over variabele.
Als casus nemen we een ficitieve populatie van 50.000 bankklanten, waarvan we de volgende variabelen weten:
sex - het geslacht van de klant
age - de leeftijd
bank - de naam van de bank
We gaan er vanuit dat dit de hele populatie is (alle klanten van twee banken binnen de onderzochte leeftijdscategorie). In werkelijkheid zal deze populatie onbekend zijn. Organisaties delen dit soort gegevens niet graag.
library(ggplot2)
N <- 50000
sex <- rep(c("M","F"), each=N/2)
age <- rep(18:22, each=N/5)
bank_m <- rep(c("ABN", "RABO"), c(N*0.3,N*0.2))
bank_f <- rep(c("ABN", "RABO"), c(N*0.1,N*0.4))
bank <- c(bank_m, bank_f)
myDF <- data.frame(sex, age, bank)
myDF[sample(1:N,10),]
summary(myDF)
sex age bank
F:25000 Min. :18 ABN :20000
M:25000 1st Qu.:19 RABO:30000
Median :20
Mean :20
3rd Qu.:21
Max. :22
tbl1 <- table(bank)
prop.tbl1 <- prop.table(tbl1)
p <- prop.tbl1[1]
tbl2 <- table(sex,bank)
tbl2
bank
sex ABN RABO
F 5000 20000
M 15000 10000
prop.table(tbl2, 1)
bank
sex ABN RABO
F 0.2 0.8
M 0.6 0.4
prop.table(tbl2, 2)
bank
sex ABN RABO
F 0.2500000 0.6666667
M 0.7500000 0.3333333
plot1 <- ggplot(myDF) +
geom_bar(aes(x=sex, fill=bank))
ggplotly(plot1)
plot2 <- ggplot(myDF) +
geom_bar(aes(x=bank, fill=sex))
ggplotly(plot2)
We bouwen nu 3 sets van steekproeven. Iedere set bevat 20 steekproeven. Per set verschilt de steekproefomvang: 10, 50 en 100. Voor iedere steekproef bepalen we de proportie ABN in het totaal. Voor iedere set steekproeven bepalen we de gemiddelde proportie.
generate_samples <- function(N = 1000, n = 10, m = 10) {
X = matrix(sample(1:N, n*m), n, m)
return(X)
}
n = 10
m = 20
aselect <- generate_samples(N = N, n = n, m = m)
samples <- sapply(1:m, function(x) myDF[aselect[,x],]$bank)
testABN <- sapply(1:m, function(x) (samples[,x] == "ABN")*1)
propABN <- sapply(1:m, function(x) sum(testABN[,x])/n)
h <- hist(propABN, xlim=range(0,1), col = "mediumaquamarine", main = paste("Histogram of m = " , m, " samples of n = ", n, "items"))
xfit<-seq(min(0),max(1),length=100)
yfit<-dnorm(xfit,mean=mean(propABN),sd=sd(propABN))
yfit <- yfit*diff(h$mids[1:2])*length(propABN)
lines(xfit, yfit, col="blue", lwd=2)

mean(propABN)
[1] 0.375
sd(propABN)
[1] 0.1585294
(p*(1-p)/n)^0.5
ABN
0.1549193
n = 50
m = 20
aselect <- generate_samples(N = N, n = n, m = m)
samples <- sapply(1:m, function(x) myDF[aselect[,x],]$bank)
testABN <- sapply(1:m, function(x) (samples[,x] == "ABN")*1)
propABN <- sapply(1:m, function(x) sum(testABN[,x])/n)
h <- hist(propABN, xlim=range(0,1), col = "tomato", main = paste("Histogram of m = " , m, " samples of n = ", n, "items"))
xfit<-seq(min(0),max(1),length=100)
yfit<-dnorm(xfit,mean=mean(propABN),sd=sd(propABN))
yfit <- yfit*diff(h$mids[1:2])*length(propABN)
lines(xfit, yfit, col="blue", lwd=2)

mean(propABN)
[1] 0.415
sd(propABN)
[1] 0.07646809
(p*(1-p)/n)^0.5
ABN
0.06928203
n = 100
m = 20
aselect <- generate_samples(N = N, n = n, m = m)
samples <- sapply(1:m, function(x) myDF[aselect[,x],]$bank)
testABN <- sapply(1:m, function(x) (samples[,x] == "ABN")*1)
propABN <- sapply(1:m, function(x) sum(testABN[,x])/n)
h <- hist(propABN, xlim=range(0,1), col = "tomato", main = paste("Histogram of m = " , m, " samples of n = ", n, "items"))
xfit<-seq(min(0),max(1),length=100)
yfit<-dnorm(xfit,mean=mean(propABN),sd=sd(propABN))
yfit <- yfit*diff(h$mids[1:2])*length(propABN)
lines(xfit, yfit, col="blue", lwd=2)

mean(propABN)
[1] 0.3905
sd(propABN)
[1] 0.04570788
(p*(1-p)/n)^0.5
ABN
0.04898979
Uit dit experiment blijkt volgende:
- Het gemiddelde (
mean) van de proportie binnen de afzonderlijke sets steekproeven ligt dicht bij de werkelijke proportie: 0.4
- Naarmate \(n\) toeneemt, neemt de spreiding rondom het gemiddelde af; er kan nauwkeuriger worden geschat wat de werkelijke waarde van het gemiddelde is.
Stel we kennen de werkelijke proportie ABN niet. De ABN zelf schat het marktaandeel op p=0.6 We kunnen op basis van de laatste set steekproeven (m=20 en n=100) onderzoeken in hoeverre de schatting realistisch is gegeven de data.
pDist <- 0.6
sdDist <- (pDist*(1-pDist)/n)^0.5
cat("De verdeling van de steekproef gemiddeldes heeft een gemiddelde van", pDist, "en een standaard deviatie van", sdDist, "\n\n")
De verdeling van de steekproef gemiddeldes heeft een gemiddelde van 0.6 en een standaard deviatie van 0.04898979
pSample <- mean(propABN)
pValue <- pnorm(pSample, mean = pDist, sd = sdDist)
Als de nulhypothese is dat het marktaandeel van de ABN 0.6 bedraagt en we nemen deze voor waar aan dan is de kans dat we een steekproef trekken met een proportie 0.3905 gelijk aan 0.0009%. Hoe overtuigd ben je na je onderzoek dat de nulhypothese (en de schatting vande ABN) juist is?
LS0tCnRpdGxlOiAiUG9wdWxhdGllcyBlbiBzdGVla3Byb2V2ZW4iCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCmBgYHtyIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0UsIGluY2x1ZGU9RkFMU0UsIHBhZ2VkLnByaW50PUZBTFNFfQpvcHRpb25zKHNjaXBlbj05OTkpCmxpYnJhcnkoZ2dwbG90MikKbGlicmFyeShwbG90bHkpCmBgYAoKV2UgYm91d2VuIGVlbiBwb3B1bGF0aWUgbWV0IGVlbiBzZXQgYWFuIHZhcmlhYmVsZW4gd2FhcnZhbiBoZXQgZ2VtaWRkZWxkZSBlbiBzcHJlaWRpbmcgYmVrZW5kIGlzLiBIaWVydWl0IGdhYW4gd2UgZGFkZWxpamsgc3RlZWtwcm9ldmVuIChzYW1wbGVzKSB0cmVra2VuIGVuIG9uZGVyem9la2VuIHdlbGsgZWZmZWN0IGRlIG9tdmFuZyB2YW4gZGUgc3RlZWtwcm9lZiBoZWVmdCBvcCBkZSBuYXV3a2V1cmlnaGVpZCB3YWFybWVlIHdlIGVlbiB1aXRzcHJhYWsga3VubmVuIGRvZW4gb3ZlciB2YXJpYWJlbGUuCgpBbHMgY2FzdXMgbmVtZW4gd2UgZWVuIGZpY2l0aWV2ZSBwb3B1bGF0aWUgdmFuIDUwLjAwMCBiYW5ra2xhbnRlbiwgd2FhcnZhbiB3ZSBkZSB2b2xnZW5kZSB2YXJpYWJlbGVuIHdldGVuOgoKLSBgc2V4YCAtIGhldCBnZXNsYWNodCB2YW4gZGUga2xhbnQKLSBgYWdlYCAtIGRlIGxlZWZ0aWpkCi0gYGJhbmtgIC0gZGUgbmFhbSB2YW4gZGUgYmFuawoKV2UgZ2FhbiBlciB2YW51aXQgZGF0IGRpdCBkZSBoZWxlIHBvcHVsYXRpZSBpcyAoYWxsZSBrbGFudGVuIHZhbiB0d2VlIGJhbmtlbiBiaW5uZW4gZGUgb25kZXJ6b2NodGUgbGVlZnRpamRzY2F0ZWdvcmllKS4gSW4gd2Vya2VsaWpraGVpZCB6YWwgZGV6ZSBwb3B1bGF0aWUgb25iZWtlbmQgemlqbi4gT3JnYW5pc2F0aWVzIGRlbGVuIGRpdCBzb29ydCBnZWdldmVucyBuaWV0IGdyYWFnLgoKYGBge3IgbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GQUxTRX0KbGlicmFyeShnZ3Bsb3QyKQpOIDwtIDUwMDAwCnNleCA8LSByZXAoYygiTSIsIkYiKSwgZWFjaD1OLzIpCmFnZSA8LSByZXAoMTg6MjIsIGVhY2g9Ti81KQpiYW5rX20gPC0gcmVwKGMoIkFCTiIsICJSQUJPIiksIGMoTiowLjMsTiowLjIpKQpiYW5rX2YgPC0gcmVwKGMoIkFCTiIsICJSQUJPIiksIGMoTiowLjEsTiowLjQpKQpiYW5rIDwtIGMoYmFua19tLCBiYW5rX2YpCm15REYgPC0gZGF0YS5mcmFtZShzZXgsIGFnZSwgYmFuaykKbXlERltzYW1wbGUoMTpOLDEwKSxdCnN1bW1hcnkobXlERikKCnRibDEgPC0gdGFibGUoYmFuaykKcHJvcC50YmwxIDwtIHByb3AudGFibGUodGJsMSkKcCA8LSBwcm9wLnRibDFbMV0KCnRibDIgPC0gdGFibGUoc2V4LGJhbmspCnRibDIKcHJvcC50YWJsZSh0YmwyLCAxKQpwcm9wLnRhYmxlKHRibDIsIDIpCgpwbG90MSA8LSBnZ3Bsb3QobXlERikgKwogIGdlb21fYmFyKGFlcyh4PXNleCwgZmlsbD1iYW5rKSkKZ2dwbG90bHkocGxvdDEpCgpwbG90MiA8LSBnZ3Bsb3QobXlERikgKwogIGdlb21fYmFyKGFlcyh4PWJhbmssIGZpbGw9c2V4KSkKZ2dwbG90bHkocGxvdDIpCmBgYAoKV2UgYm91d2VuIG51IDMgc2V0cyB2YW4gc3RlZWtwcm9ldmVuLiBJZWRlcmUgc2V0IGJldmF0IDIwIHN0ZWVrcHJvZXZlbi4gUGVyIHNldCB2ZXJzY2hpbHQgZGUgc3RlZWtwcm9lZm9tdmFuZzogMTAsIDUwIGVuIDEwMC4gVm9vciBpZWRlcmUgc3RlZWtwcm9lZiBiZXBhbGVuIHdlIGRlIHByb3BvcnRpZSBBQk4gaW4gaGV0IHRvdGFhbC4gVm9vciBpZWRlcmUgc2V0IHN0ZWVrcHJvZXZlbiBiZXBhbGVuIHdlIGRlIGdlbWlkZGVsZGUgcHJvcG9ydGllLgoKYGBge3J9CmdlbmVyYXRlX3NhbXBsZXMgPC0gZnVuY3Rpb24oTiA9IDEwMDAsIG4gPSAxMCwgbSA9IDEwKSB7CiAgWCA9IG1hdHJpeChzYW1wbGUoMTpOLCBuKm0pLCBuLCBtKQogIHJldHVybihYKQp9CgpuID0gMTAKbSA9IDIwCmFzZWxlY3QgPC0gZ2VuZXJhdGVfc2FtcGxlcyhOID0gTiwgbiA9IG4sIG0gPSBtKQoKc2FtcGxlcyA8LSBzYXBwbHkoMTptLCBmdW5jdGlvbih4KSBteURGW2FzZWxlY3RbLHhdLF0kYmFuaykKdGVzdEFCTiA8LSBzYXBwbHkoMTptLCBmdW5jdGlvbih4KSAoc2FtcGxlc1sseF0gPT0gIkFCTiIpKjEpCnByb3BBQk4gPC0gc2FwcGx5KDE6bSwgZnVuY3Rpb24oeCkgc3VtKHRlc3RBQk5bLHhdKS9uKQpoIDwtIGhpc3QocHJvcEFCTiwgeGxpbT1yYW5nZSgwLDEpLCBjb2wgPSAibWVkaXVtYXF1YW1hcmluZSIsIG1haW4gPSBwYXN0ZSgiSGlzdG9ncmFtIG9mIG0gPSAiICwgbSwgIiBzYW1wbGVzIG9mIG4gPSAiLCBuLCAiaXRlbXMiKSkKCnhmaXQ8LXNlcShtaW4oMCksbWF4KDEpLGxlbmd0aD0xMDApIAp5Zml0PC1kbm9ybSh4Zml0LG1lYW49bWVhbihwcm9wQUJOKSxzZD1zZChwcm9wQUJOKSkgCnlmaXQgPC0geWZpdCpkaWZmKGgkbWlkc1sxOjJdKSpsZW5ndGgocHJvcEFCTikgCmxpbmVzKHhmaXQsIHlmaXQsIGNvbD0iYmx1ZSIsIGx3ZD0yKQoKbWVhbihwcm9wQUJOKQpzZChwcm9wQUJOKQoocCooMS1wKS9uKV4wLjUKYGBgCgpgYGB7cn0KbiA9IDUwCm0gPSAyMAphc2VsZWN0IDwtIGdlbmVyYXRlX3NhbXBsZXMoTiA9IE4sIG4gPSBuLCBtID0gbSkKCnNhbXBsZXMgPC0gc2FwcGx5KDE6bSwgZnVuY3Rpb24oeCkgbXlERlthc2VsZWN0Wyx4XSxdJGJhbmspCnRlc3RBQk4gPC0gc2FwcGx5KDE6bSwgZnVuY3Rpb24oeCkgKHNhbXBsZXNbLHhdID09ICJBQk4iKSoxKQpwcm9wQUJOIDwtIHNhcHBseSgxOm0sIGZ1bmN0aW9uKHgpIHN1bSh0ZXN0QUJOWyx4XSkvbikKaCA8LSBoaXN0KHByb3BBQk4sIHhsaW09cmFuZ2UoMCwxKSwgY29sID0gInRvbWF0byIsIG1haW4gPSBwYXN0ZSgiSGlzdG9ncmFtIG9mIG0gPSAiICwgbSwgIiBzYW1wbGVzIG9mIG4gPSAiLCBuLCAiaXRlbXMiKSkKCnhmaXQ8LXNlcShtaW4oMCksbWF4KDEpLGxlbmd0aD0xMDApIAp5Zml0PC1kbm9ybSh4Zml0LG1lYW49bWVhbihwcm9wQUJOKSxzZD1zZChwcm9wQUJOKSkgCnlmaXQgPC0geWZpdCpkaWZmKGgkbWlkc1sxOjJdKSpsZW5ndGgocHJvcEFCTikgCmxpbmVzKHhmaXQsIHlmaXQsIGNvbD0iYmx1ZSIsIGx3ZD0yKQoKbWVhbihwcm9wQUJOKQpzZChwcm9wQUJOKQoocCooMS1wKS9uKV4wLjUKYGBgCgpgYGB7cn0KbiA9IDEwMAptID0gMjAKYXNlbGVjdCA8LSBnZW5lcmF0ZV9zYW1wbGVzKE4gPSBOLCBuID0gbiwgbSA9IG0pCgpzYW1wbGVzIDwtIHNhcHBseSgxOm0sIGZ1bmN0aW9uKHgpIG15REZbYXNlbGVjdFsseF0sXSRiYW5rKQp0ZXN0QUJOIDwtIHNhcHBseSgxOm0sIGZ1bmN0aW9uKHgpIChzYW1wbGVzWyx4XSA9PSAiQUJOIikqMSkKcHJvcEFCTiA8LSBzYXBwbHkoMTptLCBmdW5jdGlvbih4KSBzdW0odGVzdEFCTlsseF0pL24pCmggPC0gaGlzdChwcm9wQUJOLCB4bGltPXJhbmdlKDAsMSksIGNvbCA9ICJ0b21hdG8iLCBtYWluID0gcGFzdGUoIkhpc3RvZ3JhbSBvZiBtID0gIiAsIG0sICIgc2FtcGxlcyBvZiBuID0gIiwgbiwgIml0ZW1zIikpCgp4Zml0PC1zZXEobWluKDApLG1heCgxKSxsZW5ndGg9MTAwKSAKeWZpdDwtZG5vcm0oeGZpdCxtZWFuPW1lYW4ocHJvcEFCTiksc2Q9c2QocHJvcEFCTikpIAp5Zml0IDwtIHlmaXQqZGlmZihoJG1pZHNbMToyXSkqbGVuZ3RoKHByb3BBQk4pIApsaW5lcyh4Zml0LCB5Zml0LCBjb2w9ImJsdWUiLCBsd2Q9MikKCm1lYW4ocHJvcEFCTikKc2QocHJvcEFCTikKKHAqKDEtcCkvbileMC41CmBgYAoKVWl0IGRpdCBleHBlcmltZW50IGJsaWprdCB2b2xnZW5kZToKCi0gSGV0IGdlbWlkZGVsZGUgKGBtZWFuYCkgdmFuIGRlIHByb3BvcnRpZSBiaW5uZW4gZGUgYWZ6b25kZXJsaWprZSBzZXRzIHN0ZWVrcHJvZXZlbiBsaWd0IGRpY2h0IGJpaiBkZSB3ZXJrZWxpamtlIHByb3BvcnRpZTogYHIgcGAKLSBOYWFybWF0ZSAkbiQgdG9lbmVlbXQsIG5lZW10IGRlIHNwcmVpZGluZyByb25kb20gaGV0IGdlbWlkZGVsZGUgYWY7IGVyIGthbiBuYXV3a2V1cmlnZXIgd29yZGVuIGdlc2NoYXQgd2F0IGRlIHdlcmtlbGlqa2Ugd2FhcmRlIHZhbiBoZXQgZ2VtaWRkZWxkZSBpcy4KClN0ZWwgd2Uga2VubmVuIGRlIHdlcmtlbGlqa2UgcHJvcG9ydGllIGBBQk5gIG5pZXQuIERlIEFCTiB6ZWxmIHNjaGF0IGhldCBtYXJrdGFhbmRlZWwgb3AgcD0wLjYgV2Uga3VubmVuIG9wIGJhc2lzIHZhbiBkZSBsYWF0c3RlIHNldCBzdGVla3Byb2V2ZW4gKG09MjAgZW4gbj0xMDApIG9uZGVyem9la2VuIGluIGhvZXZlcnJlIGRlIHNjaGF0dGluZyByZWFsaXN0aXNjaCBpcyBnZWdldmVuIGRlIGRhdGEuCgpgYGB7cn0KcERpc3QgPC0gMC42CgpzZERpc3QgPC0gKHBEaXN0KigxLXBEaXN0KS9uKV4wLjUKY2F0KCJEZSB2ZXJkZWxpbmcgdmFuIGRlIHN0ZWVrcHJvZWYgZ2VtaWRkZWxkZXMgaGVlZnQgZWVuIGdlbWlkZGVsZGUgdmFuIiwgcERpc3QsICJlbiBlZW4gc3RhbmRhYXJkIGRldmlhdGllIHZhbiIsIHNkRGlzdCwgIlxuXG4iKQoKcFNhbXBsZSA8LSBtZWFuKHByb3BBQk4pCgpwVmFsdWUgPC0gcG5vcm0ocFNhbXBsZSwgbWVhbiA9IHBEaXN0LCBzZCA9IHNkRGlzdCkKCmBgYAoKQWxzIGRlIG51bGh5cG90aGVzZSBpcyBkYXQgaGV0IG1hcmt0YWFuZGVlbCB2YW4gZGUgQUJOIGByIHBEaXN0YCBiZWRyYWFndCBlbiB3ZSBuZW1lbiBkZXplIHZvb3Igd2FhciBhYW4gZGFuIGlzIGRlIGthbnMgZGF0IHdlIGVlbiBzdGVla3Byb2VmIHRyZWtrZW4gbWV0IGVlbiBwcm9wb3J0aWUgYHIgcFNhbXBsZWAgZ2VsaWprIGFhbiBgciByb3VuZChwVmFsdWUsIDYpKjEwMGAlLiBIb2Ugb3ZlcnR1aWdkIGJlbiBqZSBuYSBqZSBvbmRlcnpvZWsgZGF0IGRlIG51bGh5cG90aGVzZSAoZW4gZGUgc2NoYXR0aW5nIHZhbmRlIEFCTikganVpc3QgaXM/Cg==